function [X,group] = compStatistics(varargin)
% FUNCTION COMPSTATISTICS(FIELDNAME, SNTHRESHOLD,CONNECTION1, CONNECTION2,...) compares
% connection types. The parameters to be compared is indicated in
% FIELDNAME. This program first checks the normality of random variables in
% each group. If it fails to reject the normality, anova is used to compare
% the groups.
% written by Taro Kiritani, tarokiritani2008@u.northwestern.edu
% 11/07/2010

snThreshold = varargin{1,2};
numOfConnectionType = nargin - 2;
X ={};
group = {};
snRatio = {};
XStat = [];
groupStat = {};
for k = 1:numOfConnectionType
    if isvector(getfield(varargin{1,k+2}(1),varargin{1,1})) % compare vectors.
        for m = 1:length(varargin{1,k + 2})
            X{k}(m) = {getfield(varargin{1,k+2}(m),varargin{1,1})};
            group{k}(m) = {getfield(varargin{1,k+2}(m),'celltype')};
            snRatio{k}(m) = {getfield(varargin{1,k+2}(m),'snRatio')};
        end
%     else % numbers are compared.
%         for m = 1:length(varargin{1,k + 1})
%             X{k} = getfield(varargin{1,k+1}(m),varargin{1,1});
% %             group{length(X)} = getfield(varargin{1,k+1}(m),'celltype');
%         end
    end    
    
    % statistical analysis will be applied to connections with >snThreshold S/N
    % ratio.
    if length(cell2mat(X{k})) == length(group) % compare numbers
        Xsn = cell2mat(X{k}) .* (cell2mat(snRatio{k}) > snThreshold);
        snIndex = find(Xsn);
        Xsn = Xsn(snIndex);   
        group{k} = group{k}(snIndex);
        groupStat = [groupStat,group{k}];
        checkNorm(Xsn)
        XStat = [XStat,Xsn];
    else %compare vectors
        Xsn = cell2mat(X{k});
        Xsn = reshape(Xsn,length(Xsn)/m,m)';
        snMat = [];
        for n = 1:size(Xsn,2)
            snMat = [snMat,(cell2mat(snRatio{k}) > snThreshold)'];
        end
        Xsn = Xsn .* snMat;
        Xsn = Xsn(find(snMat(:,1)),:);
        XStat = [XStat;Xsn];
        group{k} = group{k}(find(snMat(:,1)));
        groupStat = [groupStat,group{k}];
        
    end
         
end


% assuming that XStat has normal distribution, use anova.
if isvector(XStat)
    [p,table, stats] = anova1(XStat,groupStat);
    c = multcompare(stats);
else
    for k = 2:size(XStat,2)
        [p,table, stats] = anova1(XStat(:,k),groupStat);
        c = multcompare(stats,'ctype','hsd');
        title(k)
    end
end

% non-parametric treatment of data using Kruskal-Wallis test.
if isvector(XStat)
    [p,table,stats] = kruskalwallis(XStat,groupStat);
    c = multcompare(stats)
else 
    for k = 2:size(XStat,2)
        [p,table, stats] = kruskalwallis(XStat(:,k),groupStat);
        c = multcompare(stats);
        title(k)
    end
end

% % pairwise ks test.
% if isvector(XStat)
%     conType = unique(groupStat);
%     for con = 1:length(conType)
%         % indices of each connection type.
%         typeInd = find(not(cellfun('isempty',strfind(groupStat,conType{con}))));
%         compval{con} = XStat(typeInd);
%     end
%     % compare all three connection types. Should be modified so that more
%     % types can be compared.
%     [h p] = kstest2(compval{1},compval{2})
%     [h p] = kstest2(compval{1},compval{3})
%     [h p] = kstest2(compval{2},compval{3})
% else
%     conType = unique(groupStat);
%     for con = 1:length(conType)
%         typeInd = find(not(cellfun('isempty',strfind(groupStat,conType{con}))));
%         compval{con} = XStat(typeInd,:);
%     end
%     
%     % compare connection types. In this case, comparisoin is done at
%     % multiple time points.
%     for timept = 1:size(XStat,2)
%         display(timept)
%         display([conType{1},' vs ',conType{2}])
%         [h p] = kstest2(compval{1}(:,timept),compval{2}(:,timept))
% 
%         display([conType{1},' vs ',conType{3}])
%         [h p] = kstest2(compval{1}(:,timept),compval{3}(:,timept))
%         
%         display([conType{2},' vs ',conType{3}])
%         [h p] = kstest2(compval{2}(:,timept),compval{3}(:,timept))
%     end
% end

end
function checkNorm (x)
% check the normality of values in X. If X is a vector, ks test is done for
% the all values in the vector X. If X is a matrix, ks test is used for
% each column of X.
if min(size(x)) > 1 % if x is a matrix. ismatrix.m is not available in matlab 2007b.
    for k = 2:size(x,2)
       [h, p, kvalue] = kstest((x(:,k) - mean(x(:,k)))/std(x(:,k)));
        if h == 1;
            display('not from normal distribution!')
        end 
    end
else
    [h, p, kvalue] = kstest((x - mean(x))/std(x));
    if h == 1;
        display('not from normal distribution!')
    end
end

end

% for k = 1:2
%     X{k} = zeros(1,length(varargin{1,k}));
%     for m = 1:length(varargin{1,k})
%         X{k}(m) = getfield(varargin{1,k}(m),varargin{1,3});
%     end
% end
% 
% [h p] = kstest2(X{1},X{2})